% Impulse Response Functions for Simple Model 
% The Forward Guidance Puzzle, Del Negro, Giannoni and Patterson
%
% Based on Giannoni's code: Nistico_model_ant_simple
% First version: October 13, 2015 

%% ========================================================================
% Calibration

% Model parameters
beta  = 0.97;
R_ss = 1/beta;
Y = 1;
delta = 0.99;
gamma = 0.15;

T = 10; 
horizon = 30;
R = zeros(horizon,1);
h = zeros(horizon,1);
w = zeros(horizon,1);
c = zeros(horizon,1);
c_j = zeros(horizon,1);
w_j = zeros(horizon,1);
h_j = zeros(horizon,1);

% Setting the Interest Rate Path 
R(1:T-1) = R_ss;
R(T) = R_ss*delta;
R(T+1:end) = R_ss;

% Calculating steady state
h_ss = Y*R_ss/(R_ss-1+gamma);
w_ss = 0; 
c_ss = (1-beta*(1-gamma))*(w_ss + h_ss);


%% Calculating Impulse Responses

% Individual H - Define this recursively
h_j(end) = h_ss;
for i = 1:horizon-1
    h_j(horizon-i) = Y + (h_j(horizon-i+1)*(1-gamma))./R(horizon-i);
end

% Individual C and W
w_j(1) = w_ss;
for i = 1:horizon
    c_j(i) = (1-beta*(1-gamma))*(w_j(i) + h_j(i)); %c_t = (1-beta(1-gamma))*(w_t + h_t)
    w_j(i+1) = (w_j(i) + Y - c_j(i))*R(i)*(1/(1-gamma)); % w_t+1 = R/1-gamma * (w_t + Y_t - c_t)
end

% Aggregate H
h = h_j ; % This is the same at the individual level amd the aggregate

% Individual C and W
w(1) = w_ss;
for i = 1:horizon
    c(i) = (1-beta*(1-gamma))*(w(i) + h(i)); %c_t = (1-beta(1-gamma))*(w_t + h_t)
    w(i+1) = (w(i) + Y - c(i))*R(i); % w_t+1 = R/1-gamma * (w_t + Y_t - c_t)
end
w_j = w_j(1:end-1);
w = w(1:end-1);


%% Figures

% Individual
subplot(3,2,1)
fig1 = plot(1:horizon, h_j./h_ss-1 );
set(fig1, {'LineWidth'}, {2});
set(fig1, {'Color'},{'b'});
ylabel('PDV of Labor Income', 'fontsize', 20);
axis([1 horizon -0.01 0.001]);
set(gca, 'YTick', [-0.01:0.002:0.001], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_ind.eps', 'psc2');

subplot(3,2,3)
fig2 = plot(1:horizon, w_j);
set(fig2, {'LineWidth'}, {2});
set(fig2, {'Color'},{'b'});
ylabel('Wealth', 'fontsize', 20);
set(gca, 'fontsize', 15);
set(gca, 'YTick', [0:0.02:0.1], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_ind.eps', 'psc2');

subplot(3,2,5)
fig3 = plot(1:horizon, c_j-1);
set(fig3, {'LineWidth'}, {2});
set(fig3, {'Color'},{'b'});
ylabel('Counsumption', 'fontsize', 20);
axis([1 horizon -0.008 0.008]);
set(gca, 'YTick', [-0.008:0.004:0.008], 'fontsize', 15);
xlabel('Individual Response', 'fontsize', 20);

%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/c_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/c_ind.eps', 'psc2');

% Aggregate
subplot(3,2,2)
fig4 = plot(1:horizon, h./h_ss -1 );
set(fig4, {'LineWidth'}, {2});
set(fig4, {'Color'},{'b'});
axis([1 horizon -0.01 0.001]);
set(gca, 'YTick', [-0.01:0.002:0.001], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_agg.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_agg.eps', 'psc2');

subplot(3,2,4)
fig5 = plot(1:horizon, w);
set(fig5, {'LineWidth'}, {2});
set(fig5, {'Color'},{'b'});
set(gca, 'YTick', [0:0.02:0.1], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_agg.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_agg.eps', 'psc2');

subplot(3,2,6)
fig6 = plot(1:horizon, c-1);
set(fig6, {'LineWidth'}, {2});
set(fig6, {'Color'},{'b'});
axis([1 horizon -0.008 0.008]);
set(gca, 'YTick', [-0.008:0.004:0.008], 'fontsize', 15);
xlabel('Aggregate Response', 'fontsize', 20);

% saveas(gcf, 'plots/irf_negative.fig');
% saveas(gcf, 'plots/irf_negative.eps', 'psc2');

w_j_keep(:,1) = w_j; 
h_j_keep(:,1) = h_j./h_ss-1;
c_j_keep(:,1) = c_j-1;
w_keep(:,1) = w;
h_keep(:,1) = h./h_ss-1;
c_keep(:,1) = c-1;

%% Alternate gamma

gamma = 0.0000001;
% Calculating steady state
h_ss = Y*R_ss/(R_ss-1+gamma);
w_ss = 0; 
c_ss = (1-beta*(1-gamma))*(w_ss + h_ss);


%% Calculating Impulse Responses

% Individual H - Define this recursively
h_j(end) = h_ss;
for i = 1:horizon-1
    h_j(horizon-i) = Y + (h_j(horizon-i+1)*(1-gamma))./R(horizon-i);
end

% Individual C and W
w_j(1) = w_ss;
for i = 1:horizon
    c_j(i) = (1-beta*(1-gamma))*(w_j(i) + h_j(i)); %c_t = (1-beta(1-gamma))*(w_t + h_t)
    w_j(i+1) = (w_j(i) + Y - c_j(i))*R(i)*(1/(1-gamma)); % w_t+1 = R/1-gamma * (w_t + Y_t - c_t)
end

% Aggregate H
h = h_j ; % This is the same at the individual level amd the aggregate

% Individual C and W
w(1) = w_ss;
for i = 1:horizon
    c(i) = (1-beta*(1-gamma))*(w(i) + h(i)); %c_t = (1-beta(1-gamma))*(w_t + h_t)
    w(i+1) = (w(i) + Y - c(i))*R(i); % w_t+1 = R/1-gamma * (w_t + Y_t - c_t)
end
w_j = w_j(1:end-1);
w = w(1:end-1);


w_j_keep(:,2) = w_j; 
h_j_keep(:,2) = h_j./h_ss-1;
c_j_keep(:,2) = c_j-1;
w_keep(:,2) = w;
h_keep(:,2) = h./h_ss-1;
c_keep(:,2) = c-1;


%% Figures
close all 
% Individual
subplot(3,2,1)
fig1 = plot(1:horizon, h_j_keep);
set(fig1, {'LineWidth'}, {2});
set(fig1, {'Color'},{'b';'r'});
ylabel('PDV of Labor Income', 'fontsize', 20);
axis([1 horizon -0.01 0.001]);
set(gca, 'YTick', [-0.01:0.002:0.001], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_ind.eps', 'psc2');

subplot(3,2,3)
fig2 = plot(1:horizon, w_j_keep);
set(fig2, {'LineWidth'}, {2});
set(fig2, {'Color'},{'b';'r'});
ylabel('Wealth', 'fontsize', 20);
set(gca, 'fontsize', 15);
set(gca, 'YTick', [0:0.02:0.1], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_ind.eps', 'psc2');

subplot(3,2,5)
fig3 = plot(1:horizon, c_j_keep);
set(fig3, {'LineWidth'}, {2});
set(fig3, {'Color'},{'b';'r'});
ylabel('Counsumption', 'fontsize', 20);
axis([1 horizon -0.008 0.008]);
set(gca, 'YTick', [-0.008:0.004:0.008], 'fontsize', 15);
xlabel('Individual Response', 'fontsize', 20);

%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/c_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/c_ind.eps', 'psc2');

% Aggregate
subplot(3,2,2)
fig4 = plot(1:horizon, h_keep);
set(fig4, {'LineWidth'}, {2});
set(fig4, {'Color'},{'b';'r'});
axis([1 horizon -0.01 0.001]);
set(gca, 'YTick', [-0.01:0.002:0.001], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_agg.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_agg.eps', 'psc2');

subplot(3,2,4)
fig5 = plot(1:horizon, w_keep);
set(fig5, {'LineWidth'}, {2});
set(fig5, {'Color'},{'b';'r'});
set(gca, 'YTick', [0:0.02:0.1], 'fontsize', 15);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_agg.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_agg.eps', 'psc2');

subplot(3,2,6)
fig6 = plot(1:horizon, c_keep);
set(fig6, {'LineWidth'}, {2});
set(fig6, {'Color'},{'b';'r'});
axis([1 horizon -0.008 0.008]);
set(gca, 'YTick', [-0.008:0.004:0.008], 'fontsize', 15);
xlabel('Aggregate Response', 'fontsize', 20);

% saveas(gcf, 'plots/irf_two_negative.fig');
% saveas(gcf, 'plots/irf_two_negative.eps', 'psc2');


%% Figures
close all 
% % Individual
% subplot(3,2,1)
% fig1 = plot(1:horizon, h_j_keep);
% set(fig1, {'LineWidth'}, {2});
% set(fig1, {'Color'},{'b';'r'});
% ylabel('PDV of Labor Income', 'fontsize', 16);
% %saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_ind.fig');
% %saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_ind.eps', 'psc2');

subplot(2,2,1)
fig2 = plot(1:horizon, w_j_keep);
set(fig2, {'LineWidth'}, {2});
set(fig2, {'Color'},{'b';'r'});
ylabel('Wealth', 'fontsize', 12);
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_ind.eps', 'psc2');

subplot(2,2,3)
fig3 = plot(1:horizon, c_j_keep);
set(fig3, {'LineWidth'}, {2});
set(fig3, {'Color'},{'b';'r'});
ylabel('Counsumption', 'fontsize', 12);
xlabel('Individual Response', 'fontsize', 15);

%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/c_ind.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/c_ind.eps', 'psc2');

% Aggregate
% subplot(3,2,2)
% fig4 = plot(1:horizon, h_keep);
% set(fig4, {'LineWidth'}, {2});
% set(fig4, {'Color'},{'b';'r'});
% %saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_agg.fig');
% %saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/h_agg.eps', 'psc2');

subplot(2,2,2)
fig5 = plot(1:horizon, w_keep);
set(fig5, {'LineWidth'}, {2});
set(fig5, {'Color'},{'b';'r'});
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_agg.fig');
%saveas(gcf, '/Users/cpatterson/Dropbox/DSGE/Long_Rate/Nistico_model/Christina_Notes/plots/w_agg.eps', 'psc2');

subplot(2,2,4)
fig6 = plot(1:horizon, c_keep);
set(fig6, {'LineWidth'}, {2});
set(fig6, {'Color'},{'b';'r'});
xlabel('Aggregate Response', 'fontsize', 15);

% saveas(gcf, 'plots/irf_two_negative.fig');
saveas(gcf, 'plots/irf_two_negative.eps', 'epsc2');
